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Abstract 

Motivated by the recent and growing interest in smart grid technology, we study the operation of DC/ AC inverters in an 
inductive microgrid. We show that a network of loads and DC/AC inverters equipped with power-frequency droop controllers 
can be cast as a Kuramoto model of phase-coupled oscillators. This novel description, together with results from the theory 
of coupled oscillators, allows us to characterize the behavior of the network of inverters and loads. Specifically, we provide a 
necessary and sufficient condition for the existence of a synchronized solution that is unique and locally exponentially stable. 
We present a selection of controller gains leading to a desirable sharing of power among the inverters, and specify the set 
of loads which can be serviced without violating given actuation constraints. Moreover, we propose a distributed integral 
controller based on averaging algorithms, which dynamically regulates the system frequency in the presence of a time- varying 
load. Remarkably, this distributed- averaging integral controller has the additional property that it preserves the power sharing 
properties of the primary droop controller. Our results hold without assumptions on identical line characteristics or voltage 
magnitudes. 

Key words: inverters; power-system control, smart power applications, synchronization, coupled oscillators, Kuramoto 
model, distributed control. 



1 Introduction 

A microgrid is a low-voltage electrical network, hetero- 
geneously composed of distributed generation, storage, 
load, and managed autonomously from the larger pri- 
mary network. Microgrids are able to connect to the 
wide area electric power system (WAEPS) through a 
Point of Common Coupling (PCC), but are also able to 
"island" themselves and operate independently. Energy 
generation within a microgrid can be highly heteroge- 
neous, any many these sources generate either variable 
frequency AC power (wind) or DC power (solar), and are 
therefore interfaced with a synchronous AC microgrid 
via power electronic devices called DC/ AC (or AC/ AC) 
power converters, or simply inverters. In islanded oper- 
ation, inverters are operated as voltage sourced inverters 
(VSIs), which act much like ideal voltage sources. It is 
through these VSIs that actions must be taken to en- 
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Fig. 1. Schematic of inverters operating in parallel. 

sure synchronization, security, power balance and load 
sharing in the network. 

Literature Review: A key topic of interest within the 
microgrid community is that of accurately sharing both 
active and reactive power among a bank of inverters op- 
erated in parallel. Such a network is depicted in Figure 
1, in which each inverter transmits power directly to 
a common load. Although several control architectures 
have been proposed to solve this problem, the so-called 
"droop" controllers have attracted the most attention, 
as they are ostensibly decentralized. The original refer- 
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Fig. 2. Mechanical analog of a Kuramoto oscillator network. 
The particles have no inertia and do not collide with another. 

ence for this methodology is [6], where Chandorkar et 
al. introduce what we will refer to as the conventional 
droop controller. For inductive lines, the droop controller 
attempts to emulate the behavior of a classical syn- 
chronous generator by imposing an inverse relation at 
each inverter between frequency and active power injec- 
tion [18]. Under other network conditions, the controller 
takes different forms [16,33,35]. Some representative ref- 
erences for the basic methodology are [30,2,21,22,20] and 
[17]. Small-signal stability analyses for two inverters op- 
erating in parallel are presented under various assump- 
tions in [10,11,23,25] and the references therein. The re- 
cent work [34] highlights some drawbacks of the con- 
ventional droop method. Distributed controllers based 
on tools from synchronous generator theory and multi- 
agent systems have also been proposed for synchroniza- 
tion and power sharing. See [26,27] for a broad overview, 
and [35,29,4,32] for various works. 

Another set of literature relevant to our investigation 
is that pertaining to synchronization of phase-coupled 
oscillators, in particular the classic and celebrated Ku- 
ramoto model. A generalization of this model considers 
n > 2 coupled oscillators, each represented by a phase 
0i G S 1 (the unit circle) and a natural frequency Vti G R. 
The system of coupled oscillators obeys the dynamics 

DiOi =n i -^ _ 1 dij sin(0i -0j), i G {1, . . . , n}, (1) 

where > is the coupling strength between the os- 
cillators i and j and Di is the time constant of the i th 
oscillator. Figure 2 shows a mechanical analog of (1), 
in which the oscillators can be visualized as a group of 
n kinematic particles, constrained to rotate around the 
unit circle. The particles rotate with preferred directions 
and speeds specified by the natural frequencies f^, and 
are connected together by elastic springs of stiffness . 
The rich dynamic behavior of the system (1) arises from 
the competition between the tendency of each oscilla- 
tor to align with its natural frequency f^, and the syn- 
chronization enforcing coupling sin(^ — 0j) with its 
neighbors. We refer to the recent surveys [1,28,12] for 
applications and theoretic results. 



The Frequency- Droop Method: The frequency- droop 
method constitutes one half of the conventional droop 
method. For inductive lines, the controller balances the 
active power demands in the network by instantaneously 
changing the frequency uji of the voltage signal at the 
i th inverter according to 

u> i = LJ*-n i (P eti -P?), (2) 

where uo* is a rated frequency, P e?i is the active electrical 
power injection at bus z, and P* is the nominal active 
power injection. The parameter > is referred to as 
the droop coefficient. 

Limitations of the Literature: Despite forming the 
foundation for the operation of parallel VSIs, the 
frequency-droop control law (2) has never been subject 
to a nonlinear analysis [34]. No conditions have been 
presented under which the controller (2) leads the net- 
work to a synchronous steady state, nor have any state- 
ments been made about the convergence rate to such a 
steady state should one exist. Stability results that are 
presented rely on linearization for the special case of 
two inverters, and sometimes come packaged with extra- 
neous assumptions [22,17]. No guarantees are given in 
terms of performance. Schemes for power sharing based 
on ideas from multi-agent systems often deal directly 
with coordinating the real and reactive power injections 
of the distributed generators, and assume implicitly that 
a low level controller is bridging the gap between the 
true network physics and the desired power injections. 
Moreover, conventional schemes for frequency restora- 
tion typically rely on a combination of local integral ac- 
tion and separation of time scales, and are generally un- 
able to maintain an appropriate sharing of power among 
the inverters (see Sections 4-5). 

Contributions: The contributions of this paper are 
four-fold. First, we begin with our key observation 
that the equations governing a microgrid under the 
frequency-droop controller can be equivalently cast as a 
generalized Kuramoto model of the form (1). We present 
a necessary and sufficient condition for the existence of 
a locally exponentially stable and unique synchronized 
solution of the closed- loop, and provide a lower bound 
on the exponential convergence rate. We also state a ro- 
bustified version of our stability condition which relaxes 
the assumption of fixed voltage magnitudes and admit- 
tances. Second, we show rigorously — and without as- 
sumptions on large output impedances or identical volt- 
age magnitudes — that if the droop coefficients are se- 
lected proportionally, then power is shared among the 
units proportionally. We provide explicit bounds on the 
set of serviceable loads. Third, we propose a distributed 
"secondary" integral controller for frequency stabiliza- 
tion. Through the use of a distributed-averaging algo- 
rithm, the proposed controller dynamically regulates the 
network frequency to a nominal value, while preserv- 
ing the proportional power sharing properties of the 
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frequency-droop controller. We show that this controller 
is locally stabilizing, without relying on the classic as- 
sumption of a time-scale separation between the droop 
and integral control loops. Fourth and finally, all results 
presented extend past the classic case of a parallel topol- 
ogy of inverters and hold for generic acyclic interconnec- 
tions of inverters and loads. 

Paper Organization: The remainder of this section 
introduces some notation and reviews some fundamen- 
tal material from algebraic graph theory, power systems 
and coupled oscillator theory. In Section 2 we motivate 
the mathematical models used throughout the rest of 
the work. In Section 3 we perform a nonlinear stability 
analysis of the frequency-droop controller. Section 4 de- 
tails results on power sharing and steady state bounds 
on power injections. In Section 5 we present and an- 
alyze our distributed- aver aging integral controller. Fi- 
nally, Section 7 concludes the paper and presents direc- 
tions for future work. 

Preliminaries and Notation: 

Sets, vectors and functions: Given a finite set V, let 
|V| denote its cardinality. Given an index set X and 
a real valued ID-array {x\, . . . , xm}? diag({a^}i G x) £ 
^|x|x|x| ^ g a ggociated diagonal matrix. We denote 
the n x n identity matrix by I n . Let l n and n be 
the n-dimensional vectors of all ones and all zeros. For 
z e R n , define = {x G R n | z T x = 0} and 
sin(z) = (sin(zi), . . . , sin(z n )) T G R n . 

Algebraic graph theory: We denote by G(V,£,A) an 
undirected and weighted graph, where V is the set of 
nodes, £ C V x V is the set of edges, and A G Rl v l x l v l is 
the adjacency matrix. If a number i G {1, . . . , \£ |} and 
an arbitrary direction is assigned to each edge {i,j} G £, 
the node-edge incidence matrix B G Rl v l x l £: l is defined 
component- wise as Bki = 1 if node k is the sink node of 
edge t and as B^i = — 1 if node k is the source node of 
edge i, with all other elements being zero. For x G R' v ', 
B T x G R'^' is the vector with components xi — Xj : with 
{ij} G £. If diag({a» J -}{ ij<7 -} G f) G R^l*^! is the diagonal 
matrix of edge weights, then the Laplacian matrix is 
given by L = Bdmg{{aij}^^ e g)B T . If the graph is 
connected, then ker(£> T ) = ker(L) = span(l|yi), and 
ker(5) = for acyclic graphs. In this case, for every x G 
that is, ^2 ieV %i = 0, there exists a unique £ G R'^' 
satisfying Kirchoff's Current Law (KCL) x = Bt; [3,9]. 
The vector x is interpreted as nodal injections, with £ 
being the associated edge flows. The Laplacian matrix 
L is positive semidefinite with eigenvalues = Ai (L) < 
^2(L) < ••• < A|v|(L). We denote the Moore-Penrose 
inverse of L by L\ and we recall from [13] the identity 
W = L^L = J| V | - |^-l| V |lj^|. 

Geometry on the n-torus: The set S 1 denotes the unit 
circle, an angle is a point G S 1 , and an arc is a con- 



nected subset of S 1 . With a slight abuse of notation, 
let \9i — 62] denote the geodesic distance between two 
angles 0i,0 2 G S 1 . The n-torus T n = 8 1 x • • • x S 1 is 
the Cartesian product of n unit circles. For 7 G [0,7r/2[ 
and a given graph G(V,£, •), let A g (t) = {0 G T' v I : 
max| i?J } G £: \0i — 0j\ < 7} be the closed set of angle ar- 
rays = . . . , 6 n ) with neighboring angles 0^ and 0j, 
{i,j} G £ no further than 7 apart. 

Synchronization: Consider the first order phase-coupled 
oscillator model (1) defined on a graph G(V, £, •). A so- 
lution : R>q — » Tl v l of (1) is said to be synchronized if 
(a) there exists a constant cj sync G R such that for each 
t > 0, 6{i) = cj sync l|v| and (b) there exists a 7 G [0, 7r/2[ 
such that 0(t) G A G (j) for each t > 0. 

^4C Power Flow: Consider a synchronous AC electrical 
network with n nodes, purely inductive admittance ma- 
trix Y G jR nxn , nodal voltage magnitudes Ei > 0, and 
nodal voltage phase angles 0^ G S 1 . The active electri- 
cal power P e ^ G R injected into the network at node 
i G {1, . . . , n} is given by [18] 

Pe,i = ^i^i I ^ I sin^i - ^ ) . (3) 



2 Problem Setup for Microgrid Analysis 

Inverter Modeling: The standard approximation in the 
microgrid literature — and the one we adopt hereafter — 
is to model an inverter as a controlled voltage source 
behind a reactance. This model is widely adopted among 
experimentalists in the microgrid field. Further modeling 
explanation can be found in [15,35,31] and the references 
therein. 

Islanded Microgrid Modeling: Figure 3 depicts an is- 
landed microgrid containing both inverters and loads. 
Such an interconnection could arise by design, or spon- 
taneously in a distribution network after an islanding 
event. An appropriate model is that of a weighted graph 
G(V,£,A) with |V| = n nodes. We consider the case of 
inductive lines, and denote by Y G jR nxn the bus ad- 
mittance matrix of the network f 1 ^ We partition the set 
of nodes as V = {Vl, V/}, corresponding to loads and in- 
verters. For {i, j} G £ , Yij is the admittance of the edge 
between nodes i and j. The output impedance of the in- 
verter can be controlled to be purely susceptive, and we 
absorb its value into the line susceptances — Im(Y^) < 0, 
{it 3} G £ . To each node i G {1, . . . , n} we assign a har- 
monic voltage signal of the form Ei(t) = Ei cos(o;*t+0i), 
where cj* > is the nominal angular frequency, E{ > 



* In some applications the inverter output impedance can 
be controlled to be highly inductive and dominate over any 
resistive effects in the network [16]. In others, the control 
law (2) is inappropriate; see [16,33,35] and Section 7. 
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Fig. 3. Schematic illustration of a microgrid, with four in- 
verters (nodes V/) supplying six loads (nodes Vl) through 
an acyclic interconnection. The dotted lines between invert- 
ers represent communication links, which will be used exclu- 
sively in Section 5. 

is the voltage amplitude, and Oi G S 1 is the voltage 
phase angle. We assume each inverter has precise mea- 
surements of its rolling time-averaged active power in- 
jection P e ,i(t) and of its frequency U0i{t), see [15] for 
details regarding this estimation. The active power in- 
jection of each inverter into the network is restricted 
to the interval [0, Pi] where Pi is the rating of inverter 
i G Vj. For the special case of a parallel interconnection 
of inverters, as in Figure 1, we will let Vl = {0} and 
V/ = {l,...,n-l}. 

3 Analysis of Frequency-Droop Control 

We now connect the frequency-droop controller (2) to 
a network of first-order phase-coupled oscillators of the 
form (1). We restrict our attention to active power flows, 
and assume the voltage magnitudes Ei are fixed at ev- 
ery bus. To begin, note that by defining Di = n" 1 and 
by writing uoi = oo* + 6^ we can equivalently write the 
frequency-droop controller (2) as 

DA = P*-Pe^ ieVi, (4) 

where P* G [0, Pi] is a selected nominal value Q Note 
that Oi is the deviation of the frequency at inverter i 
from the nominal frequency cj*. Using the active load 
flow equations (3), the droop controller (4) becomes 

DA=P:-^ n j=1 E z E J \Y ZJ \sm(O z -0 J ), i G Vj. (5) 

For a constant power loads {P*}iev L we mus ^ also satisfy 
the \Vl\ power balance equations 

= P*- V .^EiEjlYijlsmiOi-ej), ieV L . (6) 



' We make no assumptions regarding the selection of droop 
coefficients. See Section 4 for more on choice of coefficients. 



If due to failure or energy shortage an inverter i G Vi is 
unable to supply power support to the network, we for- 
mally set Di = P* =0, which reduces (5) to a load as in 
(6). If the droop-controlled system (5)-(6) reaches syn- 
chronization, then we can — without loss of generality 
— transform our coordinates to a rotating frame of ref- 
erence, where the synchronization frequency is zero and 
the study of synchronization reduces to the study of equi- 
libria. In this case, it is known that the equilibrium point 
of interest for the differential-algebraic system shares the 
same stability properties as the same equilibrium of the 
corresponding singularly perturbed system [7, Theorem 
13.1], where the constant power loads P* are replaced 
by a frequency-dependent loads P* — DiOi, for i G Vl 
and for some sufficiently small Di > 0. We can now iden- 
tify a singularly perturbed droop-controlled system of 
the form (5)— (6) with a network of Kuramoto oscillators 
described by (1) and arrive at the following insightful 
relation. 

Lemma 1 (Equivalence of Perturbed Droop- 
Controlled System and Kuramoto Model). The 

following two models are equivalent: 

(i) The singularly perturbed droop- controlled network 
(5) -(6), with frequency- dependent loads P* — Di$i 
with Di > instead of constant power loads P* G R 
(i G Vl), droop coefficients rti = 1/Di > 0, nominal 
power injections P* G R (i G Vi), nodal voltage 
phases 0i G S 1 ; nodal voltages magnitudes Ei > 0, 
and bus admittance matrix Y G jW lxn . 
(ii) The generalized Kuramoto model (1), with time con- 
stants Di > 0, natural frequencies G R ; phase 
angles 0i G S 1 and coupling weights aij > 0. 

Moreover, the parametric quantities of the two models 
are related via P* = and EiEj \Yij | = . 

In light of Lemma 1 and for notational simplicity, we 
define the matrix of time constants (inverse droop coef- 
ficients) D = diag(0|v L |, {Di}i e v T ), the vector of loads 
and nominal power injections P* = (P*, . . . , P*) T , and 
for G £ we write = EiEj\Yij\. The drop- 

controlled system (5)— (6) then reads in vector nota- 
tion as 

D0 = P*- Bdmg({a ZJ } {hj}eS )sin(B T 0), (7) 

where = (9 1 , . . . , n ) T and B G R nx l £ l is the node-edge 

incidence matrix of the underlying graph G(V, 

A natural question now arises: under what conditions 

on the power injections, network topology, admittances, 

and droop coefficients does the differential- algebraic 

closed-loop system (5)-(6) possess a stable, synchronous 

solution? 

Theorem 2 (Existence and Stability of Sync'd So- 
lution). Consider the frequency -droop controlled system 
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(5) -(6) defined on an acyclic network with node-edge in- 
cidence matrix B. Define the scaled power imbalance Cc; avg 

by ^ = (Eti ^*)/(E* ev , °i) e R > and let ^ eRlsl 

be the unique vector of edge power flows satisfying KCL, 
given implicitly by P* — uj avg Dl n = B£. The following 
two statements are equivalent: 

(i) Synchronization: There exists an arc length 7 G 
[0, 7r/2[ such that the closed-loop system (5) -(6) pos- 
sess a locally exponentially stable and unique syn- 
chronized solution t ^ Q*(t) G Ac (7) for all t > 0; 
(ii) Flow Feasibility: The power flow is feasible, i.e., 

r4||diag({a ii } {lJ}e£ )- 1 C||oo<l. (8) 

If the equivalent statements (i) and (ii) hold true, then the 
quantities T G [0, 1[ andj G [0, 7r/2[ are related uniquely 
via T = sin (7), and the following statements hold: 

a) Explicit Synchronized Solution: The synchro- 
nized solution satisfies 6*(t) = 60 + (u sync tl n ) 
(mod 2tt) for some #0 G Ac* (7), where a; sync = 
^avg; and the synchronized angular differences sat- 
isfy sin(B T 6*) = dia l g({a ij } {iJ}e e)€; 

b) Explicit Synchronization Rate: The local ex- 
ponential synchronization rate is no worse than 



ML) 



max i€ v, A 



(9) 



where L = Pdiag({a^} ^jy e g)B T is the Laplacian 
matrix of the network with weights {a>ij}{i,j}e£- 

Remark 3 (Physical Interpretation) From the 
droop controller (4), it holds that P* — uj sync Dl n G 1^ 
is the vector of steady state power injections. The power 
injections therefore satisfy the Kirchoff current law, and 
£ G R'^l is the associated vector of power flows along 
edges [9]. Physically, the parametric condition (8) there- 
fore states that the active power flow along each edge 
be feasible, i.e., less than the physical maximum = 
EiEj\Yij\. While the necessity of this condition seems 
plausible, its sufficiency is perhaps surprising. Theorem 
2 shows that equilibrium power flows are invariant under 
constant scaling of all droop coefficients, as overall scal- 
ing of D appears inversely in oj avg . Although grid stress 
varies with specific application and loading, the condition 
(8) is typically satisfied with a large margin of safety - a 
practical upper bound for 7 would be 10° . 



can consider the auxiliary system associated with (7) 
defined by 

D6 = P- Pdiag({a^} { ,, j} ^)sin(P T ^), (10) 



^av g A for 



where P; = P* for i G Vl and Pi = P i 
i G V/. Since P G 1^, system (10) has the property that 
^avg = and represents the dynamics (7) in a reference 
frame rotating at an angular frequency 0J avg . Thus, fre- 
quency synchronized solutions of (7) correspond one-to- 
one with equilibrium points of the system (10). Given 
the Laplacian matrix L = Pdiag({aij}{^j} G £)P T , (10) 
can be equivalently rewritten in the insightful form 

DO = Bdmg({a lJ } {hj}eS )-(B T L^P-sin(B T 0)) . (11) 



Here, we have made use of the facts that 
PP f = L^L = I n - ±l n l£ and P G l£. 
Since ker(P) = 0, equilibria of (11) must satisfy 
P T Ltp = B T L^Bi = sin(P T (9). We claim that 
B T L^B = diag({aij}{ i5j } G £:) _1 . To see this, define 
X = B T L^B and notice that XdidLg({aij}^jy e g)B T = 
B T LHBdmg({a ij } {i)j}e£ )B T ) = B T V [ L = B T . 
Since ker(P) = 0, it therefore holds that 
Xdia,g({aij}{ijy e £) = and the result follows. 

Hence, equilibria of (11) satisfy 



diag({a ij -} {iJ}G£ ) x f = sin(P T (9). 



(12) 



Equation (12) is uniquely solvable for f G Af)), 7 G 
[0,7r/2[, if and only if Y = max^^^/a^-) < sin(7). 
Since the right-hand side of the condition V < sin(7) 
is a concave and monotonically increasing function of 
7 G [0,7r/2[, there exists an equilibrium #* G ^g{i) f° r 
some 7 G [0, 7r/2[ if and only if the condition T < sin(7) 
is true with the strict inequality sign for 7 = ir/2. This 
leads immediately to the claimed condition T < 1. In 
this case, the explicit equilibrium angles are then ob- 
tained from the n decoupled equations (12) . See [14, The- 
orems 1, 2(G1)] for additional information. Local expo- 
nential stability of the equilibrium #* G A (7) is estab- 
lished by recalling the equivalence between the index- 
1 differential- algebraic system (11) and an associated 
reduced set of pure differential equations (see also the 
proof of (b)). In summary, the above discussion shows 
the equivalence of (i) and (ii) and statement (a). 

To show statement (b), consider the linearization of the 
dynamics (10) about the equilibrium #* G A (7) given by 



PROOF. To begin, note that if a solution t \-> 0(t) 
to the system (7) is frequency synchronized, then by 
definition there exists an oj sync G R such that 9(t) = 
^syncln for all t > 0. Summing over all equations (5)- 
(6) gives cj sync = oj avg . Without loss of generality, we 



°|Vt| 




~I\v L \ " 




~L LL 


Lli 




~A8 L ~ 


_A0!_ 




DJ 1 




Lil 


Ln 




A 9 1. 



where we have partitioned the matrix L(0*) = 
Pdiag({a^- cos(0* — 9*)}^^ e g)B T according to load 
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nodes Vl and inverter nodes V/, and defined Dj = 
diag({A}ieV/)- Since #* G A G (7), the matrix L(0*) is 
a Laplacian and thus positive semidefinite with a sim- 
ple eigenvalue at zero corresponding to rotational invari- 
ance of the dynamics under a uniform shift of all an- 
gles. It can be easily verified that the upper left block 
Lll of L(Q*) is nonsingular [13], or equivalently, 0* is 
a regular equilibrium point. Solving the set of \Vl\ al- 
gebraic equations and substituting into the dynamics 
for A0/, we obtain d(A0 7 )/dt = -DJ 1 L red (9*)A9 I , 
where L red = Lu — LjlL^Llj. The matrix L re ^(6*) G 
^|Vj|x|V/| - 1S a | go a L a pi ac i an matrix, and therefore 
shares the same properties as L(0*) [13]. Thus, it is 
the second smallest eigenvalue of DJ 1 L re d(0*) which 
bounds the convergence rate of the linearization, and 
hence the local convergence rate of the dynamics (10). 
A simple bound on ^(DJ 1 L ve &(6*)) can be obtained 
via the Courant-Fischer Theorem [24]. For x G ljvjp 
let y = D X j 2 x^ and note that x T L re a(6*)x/(x T Djx) = 



(ii) Power Injection Feasibility: 



y T DJ 1/2 L red (0* )D- iri y/{y T y). Thus, y G (DJ^l^) 1 - 



-1/2 



1/2- 



is an eigenvector of D T L re a(6*)Dj 1/z with eigen- 

— 1/2 

value /i G K. if and only if x = D T y is an eigenvector of 
DJ 1 L Te &(6*) with eigenvalue ja. For y ^ 0| Vj |, we obtain 



1/2 



A 2 (i>' / ^red(c/ J J = mm = 

2/ G(^7 1/2 iiv J ,)- y T y 

. x T L red (r)x ^ 1 . x T L red (r)x 
= mm ^— > — mm - 

xGl^! X I D I X maXieV/ ArcGlj^ , ^ ^ 

> A 2 (L red (0*)) > A 2 (L(0*)) 



max iG Vj A max iGV j A ' 

where we have made use of the spectral interlacing 
property of Schur complements [13] in the final in- 
equality. Since #* G Ag(t), the eigenvalue A 2 (£(#*)) 
can be further bounded as A 2 (L(0*)) > A 2 (L)cos(7), 
where L = BdidLg({aij}{ijy e g)B T is the Laplacian 
with weights {o>ij}{ij}e£- This fact and the identity 
cos(7) = cos(sin _1 (r)) = \/l — V 2 complete the proof. 
□ 



An analogous stability result for inverters operating in 
parallel now follows as a corollary. 

Corollary 4 (Existence and Stability of Sync'd 
Solution for Parallel Inverters). Consider a parallel 
interconnection of inverters, as depicted in Figure 1. The 
following two statements are equivalent: 



T 4 max | (if - o; avg A)/flio| < 1. (13) 
ieVi 



PROOF. For the parallel topology of Figure 1 there 
is one load fed by n — 1 inverters, and the incidence 
matrix of the graph G(V,£,A) takes the form B = 
[— l n _i 7 n _i] T . Letting P be as in the previous proof, 
we note that £ is given uniquely as £ = (B T B)~ 1 B T P. 
In this case, a set of straightforward but tedious matrix 
calculations reduce condition (8) to condition (13). □ 



Our analysis so far has been based on the assumption 
that each term = EiEj\Yij \ is a constant and known 
parameter for all {i, j} G £. In a realistic power system, 
both effective line susceptances magnitudes and voltage 
magnitudes are dynamically adjusted by additional con- 
trollers. Our analysis so far has been based on the as- 
sumption that each term = EiEj\Yij\ is a constant 
and known parameter for all {i, j} G £ . In a realistic 
power system, both effective line susceptances magni- 
tudes and voltage magnitudes are dynamically adjusted 
by additional controllers. The following result states that 
as long as these controllers can regulate the effective sus- 
ceptances and nodal voltages above prespecified lower 
bounds \Yij\ and Ej, the stability results of Theorem 2 
go through with little modification. 

Corollary 5 (Robustified Stability Condition). 

Consider the frequency- droop controlled system (5) -(6). 
Assume that the nodal voltage magnitudes satisfy Ei > 
Ei > for all i G {l,...,n} ; and that the line sus- 
ceptance magnitudes satisfy |Y^| > |Y^| > for all 

{ij} G £. For {ij} G £, ~~ 
following two statements are equivalent: 



ne = EjEj\Yjj\. The 



(i) Robust Synchronization: For all possible volt- 
age magnitudes Ei > Ei and line susceptance mag- 
nitudes \Yij\ > \Yij\, there exists an arc length 
7 G [0,7r/2[ such that the closed-loop system (5)- 
(6) possess a locally exponentially stable and unique 
synchronized solution t H> 0*(t) G Ac (7) for all 
t > 0; and 

(ii) Worst Case Flow Feasibility: The active power 
flow is feasible for the worst case voltage and line 
susceptances magnitudes, that is, 

||diag({o ! ,} { ,, j} ^)- 1 £|| 00 < 1. (14) 



(i) Synchronization: There exists an arc length 7 G 
[0,7r/2[ such that the closed-loop system (7) possess 
a locally exponentially stable and unique synchro- 
nized solution t 1— >• 9*(t) G Ac (7) for all t > 0; 



PROOF. The result follows by noting that (resp. 
Q±ij) appears exclusively in the denominator of (8) (resp. 
(14)), and that the vector f G defined in Theorem 2 
does not depend on the voltages or line susceptances. □ 
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Finally, regarding the assumption of purely inductive 
lines, we note that since the eigenvalues of a matrix are 
continuous functions of its entries, the exponential sta- 
bility property established in Theorem 2 is robust, and 
the stable synchronous solution persists in the presence 
of sufficiently small line conductances [8]. This robust- 
ness towards lossy lines is also illustrated in the simula- 
tion study of Section 6. 



4 Power Sharing and Actuation Constraints 

While Theorem 2 gives the necessary and sufficient con- 
dition for the existence of a synchronized solution to the 
closed-loop system (5)-(6), it offers no immediate guid- 
ance on how to select the control parameters P* and Di 
to satisfy the actuation constraint P e ,i £ [0, Pi]. The fol- 
lowing definition gives the proper criteria for selection. 

Definition 6 (Proportional Droop Coefficients). 

The droop coefficients are selected^ proportionally if 
P*/Di = Pf/Dj and P*/P z = P*/P for all ij G Vj. 

Theorem 7 (Power Flow Constraints and Power 
Sharing). Consider a synchronized solution of the 
frequency -droop controlled system (5) -(6), and let the 
droop coefficients be selected proportionally. Define the 
total load Pi, = J2 ieVL P* . The following two statements 
are equivalent: 

(i) Injection Constraints: < P e ,i < Pi, G Vj; 
(ii) Load Constraint: — J2jev T Pj < ^Pl < 0. 

Moreover, the inverters share the total load Pl pro- 
portionally according to their power ratings, that is, 
Pe,i/ Pi = P e ,j/Pj, for each ij G Vj. 



PROOF. From (4), the steady state active power in- 
jection at each inverter is given by P e ,i = P* — ^sync^i- 
By imposing P e ^ > for each i G V/, substituting the 
expression for Co> sy nc, and rearranging terms, we obtain, 
for each i G Vj, 




where in the final equality we have used Definition 6. 
Along with the observation that P e i > if and only if 
Pe,j > (i, j G V/), this suffices to show that < P e ,i for 
each z G Vj- if and only if Pl < 0. If we now impose for 
i G Vi that P e? ^ < Pi and again use the expression for 



cj sync along with Definition 6, a similar calculation yields 
Pei <Pi ^Pl > -— V P * = - V P 7 . 

Along with the observation that P e ^ < Pi if and only if 
Pe,j < Pj (hj £ V/), this shows that P e ,i < Pi for each 
i G Vj if and only if the total load Pl satisfies the above 
inequality. In summary, we have demonstrated two if 
and only if inequalities, which when taken together show 
the equivalence of (i) and (ii). To show the final state- 
ment, note that the fraction of the rated power capacity 
injected by the i th inverter is given by 

Pe,i _ P{ — U sync Di _ Pj ~ ^sync^j _ P e ,j 
"P. "P. " "P. "P. 

1 ^ 1 % 1 j 1 j 

for each j G Vi . This completes the proof. □ 



Power sharing results for parallel inverters supplying 
a single load follow as a corollary of Theorem 7, with 
Pl = Pq • Note that the coefficients Di must be selected 
with global knowledge. The droop method therefore re- 
quires a centralized design for power sharing despite its 
decentralized implementation. We remark that Theorem 
7 holds independently of the network voltage magnitudes 
and line admittances. 

5 Distributed PI Control in Microgrids 

As is evident from the expression for co> syn c i n Theorem 
2, the frequency-droop method typically leads to a devi- 
ation of the steady state operating frequency uj* + uj sync 
from the nominal value co>*. Again in light of Theorem 2, 
it is clear that modifying the nominal active power in- 
jection P* via the transformation P* — > P* — u sync Di 
(for i G Vj) in the controller (5) will yield zero steady 
state frequency deviation (c.f. the auxiliary system (10) 
with cD S ync = 0). Unfortunately, the information to cal- 
culate o; sync is not available locally at each inverter. As 
originally proposed in [6], after the frequency of each 
inverter has converged to cj sync , a slower, "secondary" 
control loop can be used locally at each inverter. Each 
local secondary controller slowly modifies the nominal 
power injection P* until the network frequency devia- 
tion is zero. This procedure implicitly assumes that the 
measured frequency value 0i (t) is a good approximation 
of 

^sync? and relies on a separation of time-scales be- 
tween the fast, synchronization-enforcing primary droop 
controller and the slower secondary integral controller. 
This methodology is employed in [6,19,17]. For large 
droop coefficients Di, this approach can be particularly 
slow (Theorem 2 (b)), with this slow response leading 
to an inability of the method to dynamically regulate 
the network frequency in the presence of a time- varying 
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n J = 1 

P e ,i = 5^ flij sin(^ - 0j) , iGVj 



1 



«. = ^(iy-t.) 



jeVi 



_ Pj_ 
lJ 1 D, Dj 



Fig. 4. Feedback diagram for the DAPI controller. By itself, 
the upper control loop is the droop controller (4). 



troller (15)-(16) are given by 



En 
oysin^i-^), ieV L , (17) 

Afc = A* - Pi - y n (Hj sm(0i - 9j) , i e V/ , (18) 

fciPi = P* - Pi - V" . a»j sin(6> i - 6^) 

z — '.7 = 1 



Pi_ _ Pj_ 
Di Da 



ieV T . (19) 



The following theorem (see Appendix A for the proof) 
establishes local stability of the desired equilibrium of 
(17)-(19) as well as the power sharing properties of the 
DAPI controller. 



load. Moreover, these decentralized integral controllers 
destroy the power sharing properties established by the 
primary droop controller. 



In what follows, we pursue an alternative scheme for 
frequency restoration which does not implicitly rely on 
a separation of time-scales as in [6,19,17]. Assuming 
the existence of a communication network among the 
inverters, we expand on the conventional frequency- 
droop design (2) and propose the distributed- averaging 
proportional-integral (DAPI) controller 



D i 6 i = P* - Pi -Pe,i, zeVj, (15) 



where pi G M is an auxiliary power variable and 
ki > is a gain, for each i G VjQ The matrix 
L c G ]Rl Vl l x l Vj l is the Laplacian matrix corresponding 
to a weighted, undirected and connected communication 
graph G C (V/, £ c , A c ) between the inverters, see Figure 3. 
The DAPI controller (15)-(16) is depicted in Figure 4, 
and will be shown to have the following two key proper- 
ties. First of all, the controller is able to quickly regulate 
the network frequency under large and rapid variations 
in load. Secondly, the controller accomplishes this regu- 
lation while preserving the power sharing properties of 
the primary droop controller (4). 



The closed-loop dynamics resulting from the DAPI con- 



Theorem 8 (Stability of DAPI-Controlled Net- 
work). Consider an acyclic network of droop- controlled 
inverters and loads in which the inverters can commu- 
nicate through the weighted graph G c , as described by 
the closed-loop system (17) -(19) with parameters P* G 
[0,Pi],Di > and ki > fori G Vi, and connected com- 
munication Laplacian L c Gi' Vl ' x ' Vi '. The following two 
statements are equivalent: 

(i) Stability of Droop Controller: The droop con- 
trol stability condition (8) holds; 

(ii) Stability of DAPI Controller: There exists an 
arc length 7 G [0,7r/2[ such that the system (18)- 
(19) possess a locally exponentially stable and unique 
equilibrium (#*,£>*) G A g (t) x R |Vj| . 

If the equivalent statements (i) and (ii) hold true, then 
the unique equilibrium is given as in Theorem 2 (ii), 



along with 



Da 



'avg 



for i G Vi. Moreover, if the 



droop coefficients are selected proportionally, then the 
DAPI controller preserves the proportional power sharing 
property of the primary droop controller. 

Note that Theorem 8 asserts the exponential stability of 
an equilibrium of the closed- loop (17)-(19), and hence, 
a synchronization frequency Co> syn c of zero. The network 
therefore synchronizes to the nominal frequency Co>*. 



6 Simulation Study 

We now illustrate the performance of our proposed DAPI 
controller (15)-(16) and its robustness to unmodeled 
voltage dynamics (see Corollary 5) and lossy lines in a 
simulation scenario. We consider two inverters operat- 
ing in parallel and supplying a variable load. The voltage 
magnitude at each inverter is controlled via the voltage- 
droop method 



t The presented results extend to discrete time and asyn- 
chronous communication, see [5]. 



E i = E* -mi(Qe,i-Qi) , *e{l,2}, 



(20) 
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Fig. 5. DAPI controlled closed- loop (17)-(19) for two invert- 
ers supplying a load which changes at t G {2s, 4s}. 



7 Conclusions 

We have examined the problems of synchronization, 
power sharing, and secondary control among droop- 
controlled inverters by leveraging tools from the theory 
of coupled oscillators, along with ideas from classical 
power systems and multi- agent systems. An issue not 
addressed in this work is a nonlinear analysis of reac- 
tive power sharing, as an analysis of the voltage-droop 
method (20) which yields simple and physically mean- 
ingful algebraic conditions for the existence of a solution 
is difficult to perform. Moreover, in the case of strongly 
mixed line conditions Im(Y^) ~ Re(Y^), neither of the 
control laws (2) nor (20) are appropriate. A provably 
functional control strategy for general interconnections 
and line conditions is an open and exciting problem. 

Acknowledgements 



Table 1 

Parameter values for simulation in Figure 5. 

The choice of resistances corresponds to a 
resistance/reactance ratio of one half. 



Parameter 


Symbol 






Value 


Nom. Frequency 


CJ*/27T 






60 Hz 


Nom. Voltages 


Et 




[120, 122] V 


Output/Line Indue. 


Li 




0.7 


,0.5] mH 


Output/Line Resist. 


Ri 




[0. 


14, 0.1] ft 


Inv. Ratings (P) 


p: = Pi 






[2,3] kW 


Inv. Ratings (Q) 


Qi 




[l,l]kvar 


Load (P) 


P5(t) 


PS € {- 


2.5 


, -5}kW 


Load (Q) 


Qo(t) 




-1, 


-2}kvar 


cj-Droop Coeff. 


Dt 


[4,6 




10 3 W • s 


E-Droop Coef. 


rrii 


[1.1 


] x 


1Q -3 vl 


Sec. Droop Coeff. 


hi 








Comm. Graph 


Gcomm 


Two nodes, 


one edge 


Comm. Laplacian 


L c 


(1000 Ws) • 


" 1 -l" 
-1 1 



where E* > (resp. Q* G R) is the nominal voltage 
magnitude (resp. nominal reactive power injection) at in- 
verter i G {1, 2}, rrii > Ois the voltage-droop coefficient, 
and Q e ,i G R is the reactive power injection (see [18] 
for details on reactive power). The simulation parame- 
ters are reported in Table 1. Note the effectiveness of the 
proposed DAPI controller (15)-(16) in quickly regulat- 
ing the system frequency. The spikes in local frequency 
displayed in Figure 5 (c) are due to the rapid change in 
load, and can be further suppressed by increasing the 
gains ki. This additional degree of freedom allows for a 
selection of primary droop coefficients Di much smaller 
than is typical in the literature (10 3 W-s, compared to 
roughly 10 5 W-s), allowing the power injections (Figure 
5 (b)) to respond quickly during transients. 
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A Proof of Theorem 8 

Consider the closed-loop (17)-(19) arising from the 
DAPI controller (15)-(16). We formulate our problem in 
the error coordinates Pi(t) = Pi(t) — Z^co^vg, and write 
(17)-(19) in vector notation as 



V 



} \Vl\ 

Oi 



P-P. 



} \Vl\ 



K$=Pj- P eJ - (I| Vj | + L.DJ 1 )^ 



(A.l) 
(A.2) 



where we have defined P e = BAsin(B T 0) = (P eL ,P GiI ), 
A = did,g({a ij } {iJ}eg ), V = blkdiag(/| Vl ,|,£>j), K = 
diag({/^}i G yj) and partitioned the vector of power in- 
jections by load and inverter nodes as P = (Pz,,P/) T . 
Equilibria of (A.1)-(A.2) satisfy 



= 



L Wl\ 







D~ 



L Wi\ 



J| Vj | D T + L C 








L \Vi 








DJ 1 



P-Pe 



Qi 



(A.3) 

The positive semidefinite matrix Qi has one dimensional 
kernel spanned by (0\ Vl \, D t 1\ Vi \, -l|y 7 |), while Q 2 is 
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positive definite. Note however that since P — P e G 1^, 
and Q2X = (P — P e , —Dj 1; p), it holds that Q2X ^ 
ker(Qi). Thus, (A. 3) holds if and only if x = n+ | Vj |; 

that is, p = p* = 0|Vj| and P — P e = n . Equivalently, 
from Theorem 2, the latter equation is solvable for a 
unique (modulo rotational symmetry) value #* G Ac (7) 
if and only if the parametric condition (8) holds. 



To show the final statement, note from the modified pri- 
mary controller (15) that the steady state power injec- 
tion at inverter i G V/ is given by P e i = P* — Pi(t = 
00) = P* — cj avgJ Di, which is exactly the steady state 
power injection when only the primary droop controller 
(4) is used. The result then follows from Theorem (7). 
This completes the proof of Theorem 8. □ 



To establish the local exponential stability of the equi- 
librium (0*, jJ*), we linearize the DAE (A.1)-(A.2) about 
the regular fixed point and eliminate the re- 

sulting algebraic equations, as in the proof of Theo- 
rem 2. The Jacobian J(#*,p*) of the reduced system of 
ordinary differential equations can then be factored as 
J((9*,p*) = -Z~ X X, where Z = blkdiag(I| Vj |, K) and 

£red(#*) " 

DJ 1 



X 



DJ 1 

T \Vi\ L c 



Z |V/| 



Thus, the problem of local exponential stability of 
(#*,p*) reduces to the generalized eigenvalue problem 
— X1X2V = \Zv, where A G R is an eigenvalue and 
v G R 2 I V/ I is the associated eigenvector. We will pro- 
ceed via a continuity- type argument. Consider momen- 
tarily a perturbed version of Xi, denoted by Xf, ob- 
tained by adding the matrix eI\ Vl \ to the lower-right 
block of Xl, where e > 0. Then for every e > 0, X\ 
is positive definite. Defining y = Zv, we can write 
the generalized eigenvalue problem X{X2V — —XZv as 
X 2 Z~ x y = —\{Xl)~ l y. The matrices on both left and 
right of this generalized eigenvalue problem are now sym- 
metric, with X2Z~ X = blkdiag(L re d, DJ 1 K~ 1 ) having 
a simple eigenvalue at zero corresponding to rotational 
symmetry. By applying the Courant-Fischer Theorem 
to this transformed problem, we conclude, for e > 
and modulo rotational symmetry, that all eigenvalues 
are real and negative. 



Now, consider again the unperturbed case with e = 
0. Notice that the matrix X2 is positive semidefinite 
with kernel spanned by (l| Vj | , 0| Vj |) corresponding to 
rotational symmetry, while X\ is positive semidefi- 
nite with kernel spanned by (— D/l|v J | , 1|V/|)- Since 
image(L red ((9*)) = 1^, image(X 2 ) n ker(Xi) = 
{O21 Vj I }' t na t is? X<iv is never in the kernel of X\. Thus 
we conclude that ker(XiX2) = ker(X2). Now we return 
to the original eigenvalue problem in the form X1X2V = 
— XZv. Since the eigenvalues of a matrix are continu- 
ous functions of the matrix entries, and ker(XiX2) = 
ker(X2), we conclude that the number of negative eigen- 
values does not change as e — » + , and the eigenvalues 
therefore remain real and negative. Hence, the equilib- 
rium (<9*,p*) of the DAE (A.1)-(A.2) is (again, modulo 
rotational symmetry) locally exponentially stable. 
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